subroutine task2d(time,vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,&
    ggu,ggl,di,f,v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI,l1,numbersEl,START_SP_2D_FROM)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    integer elemsCount2D,typeTask,calcWK,calcVP,timeLevel,order,START_SP_2D_FROM,start
    real*8 eps2d,eps3d,eps3dNS,radiusRED,time,vertex2D,params,ggT,ggu,ggl,di,v1,v2,v3,v4,v5,&
    ggl2,ggu2,di2,f,hFI,dt,dt1,dt0
    character*4 Ext
    dimension time(nt),vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),&
    v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),&
    di2(nodesNotTermCount2D),f(nodesNotTermCount2D),hFI(nodesCount2D,3),order(3),&
    l1(nodesBoundCount2D),numbersEl(4,nodesCount2D)

    order(1)=1
    order(2)=2
    order(3)=3

    if(START_SP_2D_FROM.eq.0) then
        !Inisial intensity
        call Hfi0(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,di,f,&
        v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI(1,order(1)),l1)
    elseif(START_SP_2D_FROM.eq.1) then
        open(500,file='.\calculations\Hfi\v2.0000',access='direct',recl=8,mode='read')
        call readFromFileDouble(500,hFI(1,order(1)),nodesCount2D)
        close(500)
    endif

    if(START_SP_2D_FROM.le.1) then
        !2-layered schema
        call Hfi1(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,di,f,&
        v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI(1,order(2)),hFI(1,order(1)),time(2)-time(1),l1)
        start = 2
    else
        start = START_SP_2D_FROM
        write(Ext,'(i4.4)') start-2
        open(500,file='.\calculations\Hfi\v2.'//Ext,access='direct',recl=8,mode='read')
        call readFromFileDouble(500,hFI(1,order(1)),nodesCount2D)
        close(500)

        write(Ext,'(i4.4)') start-1
        open(500,file='.\calculations\Hfi\v2.'//Ext,access='direct',recl=8,mode='read')
        call readFromFileDouble(500,hFI(1,order(2)),nodesCount2D)
        close(500)
    endif

    !3-layered schema
    do timeLevel=start,nt-1
        dt=time(timeLevel+1)-time(timeLevel-1)
        dt0=time(timeLevel+1)-time(timeLevel)
        dt1=time(timeLevel)-time(timeLevel-1)
        call Hfi2(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,di,f,&
        v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI(1,order(3)),hFI(1,order(1)),hFI(1,order(2)),&
        dt,dt1,dt0,timeLevel,l1)
        l=order(1)
        order(1)=order(2)
        order(2)=order(3)
        order(3)=l
    enddo

    call close_files_task2d()
    print "(A)", 'end H_fi'
end
!*************************************************************************
subroutine Hfi0(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,&
    di,f,v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI,l1)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    common /files_common/ infite_0,num_slae,log_txt,typeTask_txt
    integer elemsCount2D,typeTask,calcWK,calcVP,typeTask_txt
    real*8 vertex2D,vertex3D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,hFI
    real*8 eps2d,eps3d,eps3dNS,radiusRED
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d)
    dimension igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),&
    v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),&
    v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),di2(nodesNotTermCount2D),hFI(nodesCount2D),&
    l1(nodesBoundCount2D)

    !Assembling global matrix
    call globalForHfi0(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,igT,jgT)

    !Taking into account boundary condition
    call firstConditionForHfi(di,f,l1,vertex2D,nodesNotTermCount2D,nodesCount2D,nodesBoundCount2D,radiusRED,1)

    !Solving SLAE
    call LOS3(hFI,f,ig,jg,ggl,ggu,di,nodesNotTermCount2D,m2d,&
    maxiter,eps2d,v1,v2,v3,v4,ggl2,ggu2,di2,v5,0,iter,residual) 

    call calculateValueForNotTermNodes(hFI,igT,jgT,ggT,nodesCount2D,&
    mTerm2d,nodesTermCount2D,nodesNotTermCount2D)

    !Writing solution in file
    open(500,file='.\calculations\Hfi\v2.0000',access='direct',recl=8,status='replace',mode='write')
    call writeToFileDouble(500,hFI,nodesCount2D)
    close(500)
2   format(A,i3,A,i5,' residual =',E)
    write(*,2) 'SP: end Hfi',0,': iter=',iter,residual
    write(log_txt,2)'SP: end Hfi',0,': iter=',iter,residual
end
!*************************************************************************
subroutine Hfi1(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,di,f,&
    v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI,hFI0,dt,l1)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,m2D,mTerm2d,&
    nodesNotTermCount2D,nMaterials2d
    common /files_common/ infite_0,num_slae,log_txt,typeTask_txt
    integer elemsCount2D,typeTask,calcWK,calcVP,typeTask_txt
    real*8 vertex2D,vertex3D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,&
    ggl2,ggu2,di2,hFI,hFI0
    real*8 eps2d,eps3d,eps3dNS,radiusRED,dt
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d)
    dimension igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),&
    v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),&
    v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),di2(nodesNotTermCount2D),hFI(nodesCount2D),hFI0(nodesCount2D),&
    l1(nodesBoundCount2D)

    !Assembling global matrix
    call globalForHfi1(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,igT,jgT,hFI0,dt)

    !Taking into account boundary condition
    call firstConditionForHfi(di,f,l1,vertex2D,nodesNotTermCount2D,nodesCount2D,nodesBoundCount2D,radiusRED,2)

    !Solving SLAE
    call LOS3(hFI,f,ig,jg,ggl,ggu,di,nodesNotTermCount2D,m2d,&
    maxiter,eps2d,v1,v2,v3,v4,ggl2,ggu2,di2,v5,0,iter,residual)

    call calculateValueForNotTermNodes(hFI,igT,jgT,ggT,nodesCount2D,&
    mTerm2d,nodesTermCount2D,nodesNotTermCount2D)

    !Writing solution in file
    open(500,file='.\calculations\Hfi\v2.0001',access='direct',recl=8,status='replace',mode='write')
    call writeToFileDouble(500,hFI,nodesCount2D)
    close(500)

2   format(A,i3,A,i5,' residual =',E)
    write(*,2) 'SP: end Hfi',1,': iter=',iter,residual
    write(log_txt,2)'SP: end Hfi',1,': iter=',iter,residual
end
!*************************************************************************
subroutine Hfi2(vertex2D,nVer2D,nInCat2D,params,igT,jgT,ggT,ig,jg,ggu,ggl,&
    di,f,v1,v2,v3,v4,v5,ggl2,ggu2,di2,hFI,hFI0,hFI1,dt,dt1,dt0,timeLevel,l1)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    common /files_common/ infite_0,num_slae,log_txt,typeTask_txt
    integer elemsCount2D,typeTask,timeLevel,calcWK,calcVP,typeTask_txt
    real*8 vertex2D,vertex3D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,&
    ggl2,ggu2,di2,hFI,hFI0,hFI1
    real*8 eps2d,eps3d,eps3dNS,radiusRED,dt,dt1,dt0
    character*4 Ext
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d)
    dimension igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),&
    v1(nodesNotTermCount2D),v2(nodesNotTermCount2D),v3(nodesNotTermCount2D),v4(nodesNotTermCount2D),&
    v5(nodesNotTermCount2D),ggl2(m2d),ggu2(m2d),di2(nodesNotTermCount2D),hFI(nodesCount2D),&
    hFI0(nodesCount2D),hFI1(nodesCount2D),l1(nodesBoundCount2D)

    !Assembling global matrix
    call globalForHfi2(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,igT,jgT,hFI0,hFI1,dt,dt1,dt0)

    !Taking into account boundary condition
    call firstConditionForHfi(di,f,l1,vertex2D,nodesNotTermCount2D,nodesCount2D,nodesBoundCount2D,radiusRED,2)

    !Solving SLAE
    call LOS3(hFI,f,ig,jg,ggl,ggu,di,nodesNotTermCount2D,m2d,&
    maxiter,eps2d,v1,v2,v3,v4,ggl2,ggu2,di2,v5,0,iter,residual)

    call calculateValueForNotTermNodes(hFI,igT,jgT,ggT,nodesCount2D,&
    mTerm2d,nodesTermCount2D,nodesNotTermCount2D)

    !Writing solution in file
    write(Ext,'(i4.4)') timeLevel
    open(500,file='.\calculations\Hfi\v2.'//Ext,access='direct',recl=8,status='replace',mode='write')
    call writeToFileDouble(500,hFI,nodesCount2D)
    close(500)

2   format(A,i3,A,i5,' residual =',E)
    write(*,2) 'SP: end Hfi',timeLevel,': iter=',iter,residual
    write(log_txt,2)'SP: end Hfi',timeLevel,': iter=',iter,residual
end
!*************************************************************************
subroutine globalForHfi0(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,&
    ggT,igT,jgT)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    integer elemsCount2D,typeTask,calcWK,calcVP
    real*8 vertex2D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnRZ
    real*8 g,r1,r2,z1,z2,sigma,b,eps2d,eps3d,eps3dNS,radiusRED
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d)
    dimension g(4),b(4,4),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),&
    ig(nodesNotTermCount2D+1),jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),&
    f(nodesNotTermCount2D)

    call clearVectorDouble(g,4)
    call clear(ggl,ggu,di,f,nodesNotTermCount2D,m2d)
    do i=1,elemsCount2D
        r1=vertex2D(1,nVer2D(1,i))
        r2=vertex2D(1,nVer2D(2,i))
        z1=vertex2D(2,nVer2D(1,i))
        z2=vertex2D(2,nVer2D(3,i))
        sigma=params(nInCat2D(i))
        if(r1.eq.0.d0)r1=1d-3

        !Assembling local matrix
        call localMatrixForHfi0(b,g,r1,r2,z1,z2,sigma)
        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer2D(1,i),nodesNotTermCount2D,&
        m2d,4,igT,jgT,ggT,nodesTermCount2D,mTerm2d,1)
    enddo
end
!*************************************************************************
subroutine globalForHfi1(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,ggT,igT,jgT,hFI0,dt)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    integer elemsCount2D,typeTask,calcWK,calcVP
    real*8 vertex2D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnRZ,hFI0
    real*8 g,r1,r2,z1,z2,sigma,b,eps2d,eps3d,eps3dNS,radiusRED,f_mu,dt,hLocal
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d),hLocal(4)
    dimension g(4),b(4,4),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),hFI0(nodesCount2D)

    call clear(ggl,ggu,di,f,nodesNotTermCount2D,m2d)
    do i=1,elemsCount2D
        r1=vertex2D(1,nVer2D(1,i))
        r2=vertex2D(1,nVer2D(2,i))
        z1=vertex2D(2,nVer2D(1,i))
        z2=vertex2D(2,nVer2D(3,i))
        sigma=params(nInCat2D(i))
        f_mu=4.d0*3.1415926d0*1d-7
        if(r1.eq.0.d0)r1=1d-3

        do j=1,4
            hLocal(j)=hFI0(nVer2D(j,i))
        enddo

        !Assembling local matrix
        call localMatrixForHfi1(b,g,r1,r2,z1,z2,sigma,f_mu,dt,hLocal)
        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer2D(1,i),nodesNotTermCount2D,&
        m2d,4,igT,jgT,ggT,nodesTermCount2D,mTerm2d,1)
    enddo
end
!*************************************************************************
subroutine globalForHfi2(ggl,ggu,di,f,ig,jg,vertex2D,nVer2D,nInCat2D,params,&
    ggT,igT,jgT,hFI0,hFI1,dt,dt1,dt0)
    common /commonInfo/ eps2d,eps3d,eps3dNS,radiusRED,maxiter,typeTask,nt,nMaterials,calcWK,calcVP
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    integer elemsCount2D
    integer typeTask,calcWK,calcVP
    real*8 vertex2D,params,ggT,ggu,ggl,di,f,v1,v2,v3,v4,v5,residual,ggl2,ggu2,di2,vnRZ,hFI0,hFI1
    real*8 g,r1,r2,z1,z2,sigma,b,eps2d,eps3d,eps3dNS,radiusRED,f_mu,dt,dt1,dt0,hLocal0,hLocal1
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d)
    dimension g(4),b(4,4),igT(nodesTermCount2D+1),jgT(mTerm2d),ggT(mTerm2d),ig(nodesNotTermCount2D+1),&
    jg(m2d),ggu(m2d),ggl(m2d),di(nodesNotTermCount2D),f(nodesNotTermCount2D),hFI0(nodesCount2D),&
    hFI1(nodesCount2D),hLocal0(4),hLocal1(4)

    call clear(ggl,ggu,di,f,nodesNotTermCount2D,m2d)
    do i=1,elemsCount2D
        r1=vertex2D(1,nVer2D(1,i))
        r2=vertex2D(1,nVer2D(2,i))
        z1=vertex2D(2,nVer2D(1,i))
        z2=vertex2D(2,nVer2D(3,i))
        sigma=params(nInCat2D(i))
        f_mu=4.d0*3.1415926d0*1d-7
        if(r1.eq.0.d0)r1=1d-3

        do j=1,4
            hLocal0(j)=hFI0(nVer2D(j,i))
            hLocal1(j)=hFI1(nVer2D(j,i))
        enddo

        !Assembling local matrix
        call localMatrixForHfi2(b,g,r1,r2,z1,z2,sigma,f_mu,dt,dt1,dt0,hLocal0,hLocal1)
        !Adding local matrix to global
        call AddLocal(di,ig,jg,ggl,ggu,f,b,g,nVer2D(1,i),nodesNotTermCount2D,&
        m2d,4,igT,jgT,ggT,nodesTermCount2D,mTerm2d,1)
    enddo
end
!*************************************************************************
subroutine localStiffnessMatrixForHfi(b,r1,r2,z1,z2)
    dimension b(4,4),br(2,2),cr(2,2),bz(2,2),cz(2,2)
    real*8 b,r1,r2,z1,z2,br,cr,bz,cz

    call local_bc_r_for_H_fi(br,cr,r1,r2)
    call local_bc_z(bz,cz,z1,z2)

    do i=1,4
        do j=1,4
            mi=mu(i)
            mj=mu(j)
            ni=nu(i)
            nj=nu(j)
            b(i,j)=(br(mi,mj)*cz(ni,nj)+bz(ni,nj)*cr(mi,mj))
        enddo
    enddo
end
!*************************************************************************
subroutine local_bc_r_for_H_fi(b,c,r1,r2)
    dimension b(2,2),c(2,2),m(2,2)
    real*8 b,c,r1,r2,h,m,d

    h=r2-r1
    d=r1/h

    c(1,1)=h/6.d0*(r1*2.d0+h/2.d0*1.d0)
    c(1,2)=h/6.d0*(r1*1.d0+h/2.d0*1.d0)
    c(2,1)=h/6.d0*(r1*1.d0+h/2.d0*1.d0)
    c(2,2)=h/6.d0*(r1*2.d0+h/2.d0*3.d0)

    m(1,1)=log(1.d0+1.d0/d)*(1+d)**2-d-3.d0/2.d0
    m(1,2)=log(1.d0+1.d0/d)*(-d-d**2)+d+1.d0/2.d0
    m(2,1)=log(1.d0+1.d0/d)*(-d-d**2)+d+1.d0/2.d0
    m(2,2)=log(1.d0+1.d0/d)*d**2-d+1.d0/2.d0

    b(1,1)=(r1+h/2.d0)/h+m(1,1)-1.d0
    b(1,2)=-(r1+h/2.d0)/h+m(1,2)
    b(2,1)=-(r1+h/2.d0)/h+m(2,1)
    b(2,2)=(r1+h/2.d0)/h+m(2,2)+1.d0
end
!*************************************************************************
subroutine firstConditionForHfi(d,f,l1,vertex,nNotTerm,n,nBound,radius,numberIteration)
    dimension d(nNotTerm),f(nNotTerm),vertex(2,n),l1(nBound)
    real*8 d,f,p,vertex,radius,r,maxR,maxZ,minR,minZ,t
    common /boundaryTask2D/maxR,maxZ,minR,minZ
    p=1e+30

    if(numberIteration.eq.1)then
        do i=1,nBound
            if(vertex(2,l1(i)).eq.maxZ.and.vertex(1,l1(i)).le.radius)then
                if(vertex(1,l1(i)).eq.0.d0) then
                    r=1d-3
                else
                    r=vertex(1,l1(i))
                endif
                d(l1(i))=p
                f(l1(i))=1.d0/(2.d0*3.1415926d0*r)*p
            else
                d(l1(i))=p
                f(l1(i))=0.d0
            endif
        enddo
    else
        do i=1,nBound
            d(l1(i))=p
            f(l1(i))=0.d0
        enddo
    endif
end 
!*************************************************************************
subroutine localMassMatrixForHfi(c,r1,r2,z1,z2)
    dimension c(4,4),br(2,2),cr(2,2),bz(2,2),cz(2,2)
    real*8 c,r1,r2,z1,z2,br,cr,bz,cz

    call local_bc_r_for_H_fi(br,cr,r1,r2)
    call local_bc_z(bz,cz,z1,z2)

    do i=1,4
        do j=1,4
            mi=mu(i)
            mj=mu(j)
            ni=nu(i)
            nj=nu(j)
            c(i,j)=cr(mi,mj)*cz(ni,nj)
        enddo
    enddo
end
!*************************************************************************
subroutine localMatrixForHfi0(A,F,r1,r2,z1,z2,sigma)
    dimension A(4,4),F(4)
    real*8 A,F,r1,r2,z1,z2,sigma
    call localStiffnessMatrixForHfi(A,r1,r2,z1,z2)
    do i=1,4
        do j=1,4
            A(i,j)=A(i,j)/sigma
        enddo
    enddo
end
!*************************************************************************
subroutine localMatrixForHfi1(A,F,r1,r2,z1,z2,sigma,f_mu,dt,hPrev)
    dimension A(4,4),F(4),C(4,4),hPrev(4)
    real*8 A,F,r1,r2,z1,z2,sigma,f_mu,dt,hPrev,C

    call localStiffnessMatrixForHfi(A,r1,r2,z1,z2)
    call localMassMatrixForHfi(C,r1,r2,z1,z2)
    call mult_sq(F,C,hPrev,4)
    do i=1,4
        do j=1,4
            A(i,j)=A(i,j)/sigma+C(i,j)*f_mu/dt
        enddo
        F(i)=F(i)*f_mu/dt
    enddo
end
!*************************************************************************
subroutine localMatrixForHfi2(A,F,r1,r2,z1,z2,sigma,f_mu,dt,dt1,dt0,h0,h1)
    dimension A(4,4),F(4),B(4,4),C(4,4),gg(4),h0(4),h1(4),g(4),coef(3)
    real*8 A,F,r1,r2,z1,z2,sigma,f_mu,gg,h0,h1,g,B,C,dt,dt1,dt0,coef
    coef(1)=(dt+dt0)/(dt*dt0)
    coef(2)=(dt)/(dt1*dt0)
    coef(3)=(dt0)/(dt*dt1)
    call localStiffnessMatrixForHfi(B,r1,r2,z1,z2)
    call localMassMatrixForHfi(C,r1,r2,z1,z2)
    call mult_sq(g,C,h1,4)
    call mult_sq(gg,C,h0,4)
    do i=1,4
        do j=1,4
            A(i,j)=B(i,j)/sigma+C(i,j)*f_mu*coef(1)
        enddo
        F(i)=g(i)*f_mu*coef(2)-gg(i)*f_mu*coef(3)
    enddo
end
!*************************************************************************
subroutine close_files_task2d()
    common /files_task2d/inf2tr_dat,nvkat2d_dat,rz_dat,nvtr_dat,tsize_dat,&
    ig_dat,jg_dat,gg_dat,l1_dat,sigma_txt
    integer rz_dat,tsize_dat,gg_dat,sigma_txt
    close(inf2tr_dat)
    close(nvkat2d_dat)
    close(rz_dat)
    close(nvtr_dat)
    close(tsize_dat)
    close(ig_dat)
    close(jg_dat)
    close(gg_dat)
    close(l1_dat)
    close(sigma_txt)
end
!*************************************************************************
!*************************************************************************
!*************************************************************************